home *** CD-ROM | disk | FTP | other *** search
/ PD ROM 1 / PD ROM Volume I - Macintosh Software from BMUG (1988).iso / Programming / Programming Tools / FORTRAN Routines / RKCONST.FOR < prev    next >
Encoding:
Text File  |  1985-11-19  |  6.2 KB  |  237 lines  |  [TEXT/ttxt]
  1.  
  2. $NOFLOATCALLS
  3. $NODEBUG
  4. $STORAGE:2
  5. c*************************************************************************
  6.  
  7.     subroutine const(ist,imeth,idef)
  8.  
  9. c****    subroutine supplies constants used in the runge-kutta
  10. c****    method,in the format required by RKP.
  11.  
  12.     implicit double precision (a-h,p)
  13.     dimension al(12),b(12,12)
  14.     common /coeffs/al,b,a1,a3,a4,a5,a6,a7,a8,a9,a10,a11,
  15.      &               a12,a13,er1,er3,er4,er5,er6,er7,er8,
  16.      &               er9,er10,er11,er12,er13,inum
  17.     common /cons/hmin,hmaxl,hfact,hlim,power
  18.  
  19.     do 10 j=1,12
  20.        do 5 i=1,12
  21.           b(i,j)=0.d0
  22. 5       continue
  23.        al(j)=0.d0
  24. 10     continue
  25.  
  26.        if(imeth.EQ.1) then
  27. c****  Fourth order runge-kutta method specified by
  28. c****  Felberg, E., COMPUTING,    6(1970)p61-71
  29. c      write(*,*) 'Fourth order Runge-Kutta-Fehlberg method'
  30.        al(1)=1.d0/4.d0
  31.        al(2)=3.d0/8.d0
  32.        al(3)=12.d0/13.d0
  33.        al(4)=1.d0
  34.        al(5)=1.d0/2.d0
  35.  
  36.        b(1,1)=1.d0/4.d0
  37.        b(1,2)=3.d0/32.d0
  38.        b(2,2)=9.d0/32.d0
  39.        b(1,3)=1932.d0/2197.d0
  40.        b(2,3)=-7200.d0/2197.d0
  41.        b(3,3)=7296.d0/2197.d0
  42.        b(1,4)=439.d0/216.d0
  43.        b(2,4)=-8.d0
  44.        b(3,4)=3680.d0/513.d0
  45.        b(4,4)=-845.d0/4104.d0
  46.        b(1,5)=-8.d0/27.d0
  47.        b(2,5)=2.d0
  48.        b(3,5)=-3544.d0/2565.d0
  49.        b(4,5)=1859.d0/4104.d0
  50.        b(5,5)=-11.d0/40.d0
  51.  
  52.        er1=1.d0/360.d0
  53.        er3=-384.d0/12825.d0
  54.        er4=-41743.d0/1429560.d0
  55.        er5=1.d0/50.d0
  56.        er6=2.d0/55.d0
  57.  
  58.        a1=16.d0/135.d0
  59.        a2=0.d0
  60.        a3=6656.d0/12825.d0
  61.        a4=28561.d0/56430.d0
  62.        a5=-9.d0/50.d0
  63.        a6=2.d0/55.d0
  64.  
  65.        inum=6
  66.        ist=4
  67.        if(idef.eq.0) hmaxl=0.5d0
  68.        hlim1=3.d0
  69.        end if
  70.  
  71.        if(imeth.EQ.2) then
  72. c****  Fifth order runge-kutta method specified by
  73. c****  Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 5)
  74. c      write(*,*) 'Fifth order Runge-Kutta-Verner method'
  75.        al(1)=1.d0/18.d0
  76.        al(2)=1.d0/6.d0
  77.        al(3)=2.d0/9.d0
  78.        al(4)=2.d0/3.d0
  79.        al(5)=1.d0
  80.        al(6)=8.d0/9.d0
  81.        al(7)=1.d0
  82.  
  83.        b(1,1)=1.d0/18.d0
  84.        b(1,2)=-1.d0/12.d0
  85.        b(2,2)=1.d0/4.d0
  86.        b(1,3)=-2.d0/81.d0
  87.        b(2,3)=4.d0/27.d0
  88.        b(3,3)=8.d0/81.d0
  89.        b(1,4)=40.d0/33.d0
  90.        b(2,4)=-4.d0/11.d0
  91.        b(3,4)=-56.d0/11.d0
  92.        b(4,4)=54.d0/11.d0
  93.        b(1,5)=-369.d0/73.d0
  94.        b(2,5)=72.d0/73.d0
  95.        b(3,5)=5380.d0/219.d0
  96.        b(4,5)=-12285.d0/584.d0
  97.        b(5,5)=2695.d0/1752.d0
  98.        b(1,6)=-8716.d0/891.d0
  99.        b(2,6)=656.d0/297.d0
  100.        b(3,6)=39520.d0/891.d0
  101.        b(4,6)=-416.d0/11.d0
  102.        b(5,6)=52.d0/27.d0
  103.        b(1,7)=3015.d0/256.d0
  104.        b(2,7)=-9.d0/4.d0
  105.        b(3,7)=-4219.d0/78.d0
  106.        b(4,7)=5985.d0/128.d0
  107.        b(5,7)=-539.d0/384.d0
  108.        b(7,7)=693.d0/3328.d0
  109.  
  110.        er1=33.d0/640.d0
  111.        er3=-132.d0/325.d0
  112.        er4=891.d0/2240.d0
  113.        er5=-33.d0/320.d0
  114.        er6=-73.d0/700.d0
  115.        er7=891.d0/8320.d0
  116.        er8=2.d0/35.d0
  117.  
  118.        a1=57.d0/640.d0
  119.        a3=-16.d0/65.d0
  120.        a4=1377.d0/2240.d0
  121.        a5=121.d0/320.d0
  122.        a7=891.d0/8320.d0
  123.        a8=2.d0/35.d0
  124.  
  125.        inum=8
  126.        ist=5
  127.        if(idef.eq.0) hmaxl=1.0d0
  128.        hlim1=4.d0
  129.        end if
  130.  
  131.        if(imeth.EQ.3) then
  132. c****  Seventh order runge-kutta method specified by
  133. c****  Verner J.H., SIAM J. Numer. Anal. V15,(1978),p772.(table 7)
  134. c      write(*,*) 'Seventh order Runge-Kutta-Verner method'
  135.        al(1)=1.d0/4.d0
  136.        al(2)=1.d0/12.d0
  137.        al(3)=1.d0/8.d0
  138.        al(4)=2.d0/5.d0
  139.        al(5)=1.d0/2.d0
  140.        al(6)=6.d0/7.d0
  141.        al(7)=1.d0/7.d0
  142.        al(8)=2.d0/3.d0
  143.        al(9)=2.d0/7.d0
  144.        al(10)=1.d0
  145.        al(11)=1.d0/3.d0
  146.        al(12)=1.d0
  147.  
  148.        b(1,1)=1.d0/4.d0
  149.        b(1,2)=5.d0/72.d0
  150.        b(2,2)=1.d0/72.d0
  151.        b(1,3)=1.d0/32.d0
  152.        b(3,3)=3.d0/32.d0
  153.        b(1,4)=106.d0/125.d0
  154.        b(3,4)=-408.d0/125.d0
  155.        b(4,4)=352.d0/125.d0
  156.        b(1,5)=1.d0/48.d0
  157.        b(4,5)=8.d0/33.d0
  158.        b(5,5)=125.d0/528.d0
  159.        b(1,6)=-1263.d0/2401.d0
  160.        b(4,6)=39936d0/26411d0
  161.        b(5,6)=-64125.d0/26411d0
  162.        b(6,6)=5520.d0/2401.d0
  163.        b(1,7)=37.d0/392.d0
  164.        b(5,7)=1625.d0/9408.d0
  165.        b(6,7)=-2.d0/15.d0
  166.        b(7,7)=61.d0/6720.d0
  167.        b(1,8)=17176.d0/25515.d0
  168.        b(4,8)=-47104.d0/25515.d0
  169.        b(5,8)=1325.d0/504.d0
  170.        b(6,8)=-41792.d0/25515.d0
  171.        b(7,8)=20237.d0/145800.d0
  172.        b(8,8)=4312.d0/6075.d0
  173.        b(1,9)=-23834.d0/180075.d0
  174.        b(4,9)=-77824.d0/1980825.d0
  175.        b(5,9)=-636635.d0/633864.d0
  176.        b(6,9)=254048.d0/300125.d0
  177.        b(7,9)=-183.d0/7000.d0
  178.        b(8,9)=8.d0/11.d0
  179.        b(9,9)=-324.d0/3773.d0
  180.        b(1,10)=12733.d0/7600.d0
  181.        b(4,10)=-20032.d0/5225.d0
  182.        b(5,10)=456485.d0/80256.d0
  183.        b(6,10)=-42599.d0/7125.d0
  184.        b(7,10)=339227.d0/912000d0
  185.        b(8,10)=-1029.d0/4180.d0
  186.        b(9,10)=1701.d0/1408.d0
  187.        b(10,10)=5145.d0/2432.d0
  188.        b(1,11)=-27061.d0/204120.d0
  189.        b(4,11)=40448.d0/280665.d0
  190.        b(5,11)=-1353775.d0/1197504.d0
  191.        b(6,11)=17662.d0/25515.d0
  192.        b(7,11)=-71687.d0/1166400d0
  193.        b(8,11)=98.d0/225.d0
  194.        b(9,11)=1.d0/16.d0
  195.        b(10,11)=3773.d0/11664.d0
  196.        b(1,12)=11203.d0/8680.d0
  197.        b(4,12)=-38144.d0/11935.d0
  198.        b(5,12)=2354425.d0/458304.d0
  199.        b(6,12)=-84046.d0/16275.d0
  200.        b(7,12)=673309.d0/1636800.d0
  201.        b(8,12)=4704.d0/8525.d0
  202.        b(9,12)=9477.d0/10912.d0
  203.        b(10,12)=-1029.d0/992.d0
  204.        b(12,12)=729.d0/341.d0
  205.  
  206.        er1=-1.d0/480.d0
  207.        er6=-16.d0/375.d0
  208.        er7=-2401.d0/528000.d0
  209.        er8=2401.d0/132000.d0
  210.        er9=243.d0/14080.d0
  211.        er10=-2401.d0/19200.d0
  212.        er11=-19.d0/450.d0
  213.        er12=243.d0/1760.d0
  214.        er13=31.d0/720.d0
  215.  
  216.        a1=31.d0/720.d0
  217.        a6=16.d0/75.d0
  218.        a7=16807.d0/79200.d0
  219.        a8=16807.d0/79200.d0
  220.        a9=243.d0/1760.d0
  221.        a12=243.d0/1760.d0
  222.        a13=31.d0/720.d0
  223.  
  224.        inum=13
  225.        ist=7
  226.        if(idef.eq.0) hmaxl=2.d0
  227.        hlim1=5.d0
  228.        end if
  229.  
  230.        power=1.d0/dble(ist+1)
  231.        hfact=0.5d0**(1.d0/dble(ist-1))
  232.        hlim=hlim1/hfact
  233.        if(idef.eq.0) hmin=1.d-06
  234. c      write(*,*)' Constants evaluated'
  235.        return
  236.        end
  237.